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Abstract 

This paper presents a methodology to determine the parameters to be used in 
the constitutive equations of Cohesive Zone Models employed in the simulation 
of delanrination in composite materials by means of decohesion finite elements. A 
closed-form expression is developed to define the stiffness of the cohesive layer. A 
novel procedure that allows the use of coarser meshes of decohesion elements in 
large-scale computations is also proposed. The procedure ensures that the energy 
dissipated by the fracture process is computed correctly. It is shown that coarse- 
rneshed models defined using the approach proposed here yield the same results as 
the models with finer meshes normally used for the simulation of fracture processes. 
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1 Introduction 


Delamination, or interfacial cracking between composite layers, is one of the 
most common types of damage in laminated fibre-reinforced composites due 
to their relatively weak interlaminar strengths. Delamination may arise un- 
der various circumstances, such as low velocity impacts, or bearing loads in 
structural joints. The delamination failure mode is particularly important for 
the structural integrity of composite structures because it is difficult to detect 
during inspection. 

The simulation of delamination using the finite element method (FEM) is 
normally performed by means of the Virtual Crack Closure Technique (VCCT) 
[1], or using decohesion finite elements [2]- [10] . 

The VCCT is based on the assumption that the energy released during de- 
lamination propagation equals the work required to close the crack back to its 
original position. Based on this assumption, the components of the energy re- 
lease rate are computed from the nodal forces and nodal relative displacements 
[1] . Delamination growth is predicted when a combination of the components 
of the energy release rate is equal to a critical value. 

There are some difficulties when using the VCCT in the simulation of pro- 
gressive delamination. The calculation of fracture parameters requires nodal 
variables and topological information from the nodes ahead and behind the 
crack front. Such calculations are tedious to perform and may require remesh- 
ing for crack propagation. 

The use of decohesion finite elements can overcome some of the above difficul- 
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ties. Decohesion finite elements can predict both the onset and the non-self- 
similar propagation of delamination. However, the simulation of progressive 
delamination using decohesion elements poses numerical difficulties related 
with the proper definition of the stiffness of the cohesive layer, the require- 
ment of extremely refined meshes, and the convergence difficulties associated 
with problems involving softening constitutive models. 

This work addresses the proper definition of interface stiffness. A closed-form 
expression for the interface stiffness is developed, replacing the empirical defi- 
nitions of the stiffness of the cohesive layer that are normally used. In addition, 
a solution that allows the use of coarse meshes in the simulation of delamina- 
tion is proposed. 

The kinematics and constitutive model of a previously proposed decohesion 
element [11]- [12], formulated for three-dimensional (3D) elements, is adapted 
to be used with two-dimensional (2D) finite element models. Finally, several 
simulations of specimens with and without initial cracks were performed, in 
order to demonstrate that the methodology proposed can accurately predict 
both crack initiation and propagation using coarse meshes. 


2 Selection of Cohesive Zone Model parameters 

2.1 Cohesive Zone Model approach 

Decohesion finite elements are used to model material discontinuities. Consider 
a domain D containing a material discontinuity, T^, which divides the domain 
D into two parts, D+ and D-, as shown in Figure 1. 
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[Figure 1 about here] 


Prescribed tractions, t tl are imposed on the boundary V p, whereas prescribed 
displacements are imposed on the boundary T u . The stress held inside the 
domain, a Vj , is related to the external loading and to the closing tractions 
in the material discontinuity through the equilibrium equations: 


(Tijj = 0 in Q (1) 

(TijUj = ti on T p (2) 

rrijiij = T+ = -r" = cr ij n j on T d (3) 

The formulation of the decohesion hnite elements is based in the Cohesive 
Zone Model (CZM) approach. The CZM approach [13]— [15] is one of the most 
commonly used tools to investigate interfacial fracture. The CZM approach 
assumes that a cohesive damage zone, or softening plasticity, develops near 
the crack tip. The CZM links the microstructural failure mechanism to the 
continuum fields governing bulk deformations. Thus, a CZM is characterized 
by the properties of the bulk material, the crack initiation condition, and the 
crack evolution function. 

Cohesive damage zone models relate traction, r, to displacement jumps, A, at 
an interface where a crack may occur. Damage initiation is related to the in- 
terfacial strength, r°, i.e., the maximum traction on the traction-displacement 
jump relation. When the area under the traction-displacement jump relation 
is equal to the fracture toughness, G c , the traction is reduced to zero, and 
new crack surfaces are formed. There are several types of constitutive equa- 
tions used in decohesion elements: Tvergaard and Hutchinson [16] proposed a 
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trapezoidal law, Cui an Wisnom [17] a perfectly plastic rule, Needleman first 
proposed a polynomial law, [18], and later an exponential law [19]. Goyal et al. 
[10] adopted Needleman’s exponential law to account for load reversal with- 
out restoration of the cohesive state. The law used in this paper is a bilinear 
relation between the tractions and the displacement jumps [6], [9], [11], [12], [20], 
(see Figure 2). 


[Figure 2 about here] 


2.2 Cohesive zone model and FEM 

In a finite element model using the CZM approach, the complete material 
description is separated into fracture properties captured by the constitutive 
model of the cohesive surface and the properties of the bulk material, captured 
by the continuum regions. 

There are two conditions to obtain a successful FEM simulation using CZM 
[21]: (a) The cohesive contribution to the global compliance before crack prop- 
agation should be small enough to avoid the introduction of a fictitious compli- 
ance to the model, and (b) the element size needs to be less than the cohesive 
zone length. 

(a) Stiffness of the cohesive zone model 

Different guidelines have been proposed for selecting the stiffness of the in- 
terface. Daudeville et al. [23] calculated the stiffness in terms of the thickness 
and the elastic modulus of the interface. The resin-rich interface between plies 
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is of the order of 10 _5 m. Therefore, the interface stiffness obtained from the 
Daudeville equations is very high. Zou et al. [24], based on their own expe- 
rience, proposed a value for the interface stiffness between 10 4 and 10 7 times 
the value of the interfacial strength per unit length. Camanho and Davila [9] 
obtained accurate predictions for Graphite/Epoxy specimens using a value of 
10 6 N/mm 3 . 

The effective clastic properties of the whole laminate depend on the properties 
of both the cohesive surfaces and the bulk constitutive relations of the plies. 
Although the compliance of the cohesive surfaces can contribute to the global 
deformation, its only purpose is to simulate fracture. Moreover, the elastic 
properties of the cohesive surfaces are mesh-dependent because the surface 
relations exhibit an inherent length scale that is absent in homogeneous de- 
formations [22], 

If the cohesive contribution to the compliance is not small enough compared 
to that of the volumetric constitutive relation, a stiff connection between two 
neighboring layers before delamination initiation is not assured. The effect of 
compliance of the interface on the bulk properties of a laminate is illustrated 
in the one- dimensional model shown in Figure 3. The traction continuity con- 
dition requires: 


a = E 3 e = KA (4) 

where a is the traction on the surface, t is the thickness of an adjacent sub- 
laminate, e = y is the transverse strain, K is the interface stiffness that 
relates the resulting tractions at the interface with the opening displacement 
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A, and E :i is the through-the-thickness Young’s modulus of the material. For 


a transversely isotropic material E 3 = E 2 . 


[Figure 3 about here] 


The effective strain of the composite is: 

St A A . . 

£ ‘’~j + t~ £+ j (5) 

Since the traction continuity condition requires that o = the equivalent 

Young’s modulus E e g can be written as a function of the Young’s modulus of 
the material, the mesh size, and the interface stiffness. Using equations (4) 
and (5), the effective Young’s modulus can be written as: 


Enff — E 3 


1 + El 
. 1 ^ Kt - 


(6) 


The effective elastic properties of the composite will not be affected by the 
cohesive surface whenever E 3 <C Kt, i.e: 


K = 


aE 3 

t 


( 7 ) 


where a is a parameter much larger than 1 (a > 1). However, large values of 
the interface stiffness may cause numerical problems, such as spurious oscilla- 
tions of the tractions [4], Thus, the interface stiffness should be large enough 
to provide a reasonable stiffness but small enough to avoid numerical problems 
such as spurious oscillations of the tractions in an element. 
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The ratio between the value of the Young modulus obtained with equation (6) 
and the Young modulus of the material, as a function of the parameter a is 
shown in Figure 4. For values of a greater than 50, the loss of stiffness due to 
the presence of the interface is less than 2%. 

[Figure 4 about here] 


The use of equation (7) is preferable to the guidelines presented in previous 
work [9], [23]- [24] because it results from mechanical considerations, and it pro- 
vides a sufficient stiffness (a times the stiffness of the material) while avoiding 
spurious oscillations caused by an excessively stiff interface. The values of the 
interface stiffness obtained with equation (7) and those used by other authors 
for a carbon/epoxy composite are shown in Table 1. The material’s transverse 
modulus is llGPa, its nominal interfacial strength is r° = 45MPa, and a is 
taken as 50. 


[Table 1 about here] 

Equation (7) gives a range of the interface stiffness between 10 5 and 5xl0 6 N/mm 3 
for a sub-laminate thickness between 0.125mm and 5mm. These values are 
close to the interface stiffness proposed by Camanho and Davila and the val- 
ues obtained from Zou’s guidelines (between 4.5xl0 5 and 4.5xl0 8 N/mm 3 ). 

(b) Length of the cohesive zone 

Under single-mode loading, interface damage initiates at a point where the 
traction reaches the maximum nominal interfacial strength, r°. For mixed- 
mode loading, interface damage onset is predicted by a criterion established 



in terms of the normal and shear tractions. Crack propagation occurs when the 
energy release rate reaches a critical value G c . The CZM approach prescribes 
the interfacial normal and shear tractions that resist separation and relative 
sliding at an interface. The tractions, integrated to complete separation, yield 
the fracture energy release rate, G c . The length of the cohesive zone l cz is 
defined as the distance from the crack tip to the point where the maximum 
cohesive traction is attained (see Figure 5). 

[Figure 5 about here] 


Different models have been proposed in the literature to estimate the length 
of the cohesive zone. Irwin [25] estimated the size of the plastic zone ahead 
of a crack in a ductile solid by considering the crack tip zone within which 
the von Mises equivalent stress exceeds the tensile yield stress. Dugdale [13] 
estimated the size of the yield zone ahead of a mode I crack in a thin plate of an 
elastic-perfectly plastic solid by idealizing the plastic region as a narrow strip 
extending ahead of the crack tip that is loaded by the yield traction. Barenblatt 
[14] provided an analogue for ideally brittle materials of the Dugdale plastic 
yield zone analysis. Hui [26] estimated the length of the cohesive zone for soft 
clastic solids, and Falk [21] and Rice [27] estimated the length of the cohesive 
zone as a function of the crack growth velocity. The expressions that result 
from these models for the case of plane stress are presented in Table 2. It is 
assumed that the relation between the critical stress intensity factor K c and 
the critical energy release rate G c can be expressed as Jy? = G C E. 

All of the models described above to predict the cohesive zone length l cz have 
the form: 
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G 

l cz = ME — % 

( r 0)2 


( 8 ) 


where E is the Young modulus of the material, G c is the critical energy release 
rate, r° is the maximum interfacial strength, and M is a parameter that 
depends on each model. The most commonly used models in the literature are 
Hillerborg’s model [15] and Rice’s model [27]. In these models, the parameter 
M is either close or exactly equal to unity. A summary of the different models 
commonly used in the literature, and the equivalent parameter M for plane 
stress are shown in Table 2. In this paper, Hillerborg’s model is used in the 
following analysis. 

For the case of orthotropic materials with transverse isotropy, the value of the 
Young’s modulus in equation (8) is the transverse modulus of the material, 

e 2 . 


[Table 2 about here] 


In order to obtain accurate FEM results using CZM, the tractions in the 
cohesive zone must be represented properly by the finite element spatial dis- 
cretization. The number of elements in the cohesive zone is: 

N e = l f (9) 

where l e is the mesh size in the direction of crack propagation. 

When the cohesive zone is discretized by too few elements, the fracture energy 
is not represented accurately and the model does not capture the continuum 
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field of a cohesive crack. Therefore, a minimum number of elements, N e , is 
needed in the cohesive zone to get successful FEM results. 

However, the minimum number of elements needed in the cohesive zone is not 
well established: Moes and Belytschko [28], based on the work of Carpinteri 
and Cornetti [29], suggest using more than 10 elements. However, Falk et 
al. [21] used between 2 and 5 elements in their simulations. In the parametric 
study by Davila and Camanho [30] , the minimum element length for predicting 
the delamination in a double cantilever beam (DCB) specimen was 1mm, 
which leads, using equation (8) with M = 1, to a length of the cohesive zone 
of 3.28mm. Therefore, 3 elements in the cohesive zone were sufficient to predict 
the propagation of delamination in Mode I. 


2.3 Guidelines for the selection of the parameters of the interface with coarser 
meshes 

One of the drawbacks in the use of cohesive zone models is that very foie 
meshes are needed to assure a reasonable number of elements in the cohesive 
zone. The length of the cohesive zone given by equation (8) is proportional 
to the fracture energy release rate ( G c ) and to the inverse of the square of 
the interfacial strength r°. For typical graphite-epoxy or glass-epoxy com- 
posite materials, the length of the cohesive zone is smaller than one or two 
millimeters. Therefore, according to equation (9), the mesh size required in 
order to have more than two elements in the cohesive zone should be smaller 
than half a millimeter. The computational requirements needed to analyze 
a large structure with these mesh sizes may render most practical problems 
intractable. 
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Alfano and Crisfield [8] observed that variations of the maximum interfacial 
strength do not have a strong influence in the predicted results, but that low- 
ering the interfacial strength can improve the convergence rate of the solution. 
The result of using a lower interfacial strength is that the length of the cohe- 
sive zone and the number of elements in the cohesive zone increase. Therefore, 
the representation of the softening response of the fracture process ahead of 
the crack tip is more accurate with a lower interface strength although the 
stress distribution in the regions near the crack tip might be altered [8]. 

It is possible to develop a strategy to adapt the length of the cohesive zone 
to a given mesh size. The procedure consists of determining the value r° of 
the interfacial strength required for a desired number of elements (7V f °) in the 
cohesive zone. From equations (8) and (9), the required interface strength is: 



Finally, the interfacial strength is chosen as: 


T = min {r°, r 0 } 


( 10 ) 


( 11 ) 


The interfacial strength is computed for each loading mode, replacing the value 
of the fracture toughness, G c , in equation (10) by the value corresponding to 
the loading mode. 

The effect of a reduction of the interfacial strength is to enlarge the cohesive 
zone, and thus, the model is better suited to capture the softening behaviour 
ahead of the crack tip. Moreover, if equation (7) is used to compute the in- 
terface stiffness, the interface stiffness will be large enough to assure a stiff 
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connection between the two neighboring layers and small enough to avoid 
spurious oscillations. 


3 FEM implementation 

The cohesive zone model used here is based on the model that was previously 
proposed by the authors [11]- [12]. The model has been adapted to be used in 
2D finite element models. A brief description of the model, with special focus 
on the kinematics and constitutive damage model, is presented. 

3. 1 Kinematics and constitutive relation of cohesive zone models 

The displacement jump across the interface of the material discontinuity re- 
quired in the constitutive model, [?!,] , can be obtained as a function of the 
displacement of the points located on the top and bottom side of the interface, 
uf and u ~ , respectively: 


IM = u t -Ui (12) 

where uf are the displacements with respect to the fixed Cartesian coordinate 
system. A co-rotational formulation is used in order to express the compo- 
nents of the displacement jumps with respect to the deformed interface. The 
coordinates Xi of the deformed interface can be written as [31]: 


Xi 




(13) 
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where X t are the coordinates of the undeformed interface. 


The components of the displacement jump tensor in the local coordinate sys- 
tem on the deformed interface, A m , are expressed in terms of the displacement 
field in global coordinates: 

A m = [Mil ( 14 ) 

where Q mi is the rotation tensor. 

The constitutive operator of the interface, D Jt , relates the element tractions, 
Tj, to the displacement jumps, Ap 


tj = DjiAi ( 15 ) 

A requirement of the constitutive relations of cohesive zone models is that the 
energy dissipated in the process of fracture needs to be computed accurately. 

Under single-mode loading, controlled energy dissipation is achieved assur- 
ing that the area under the traction- displacement jump relation equals the 
corresponding fracture toughness, as illustrated in Figure 6. 

[Figure 6 about here] 


Under mixed-mode loading, a criterion established in terms of an interaction 
between components of the energy release rate needs to be used to predict 
crack propagation. 


14 



The constitutive damage model used here, formulated in the context of Dam- 
age Mechanics (DM), was previously proposed by the authors [11], [12]. All 
the details of the constitutive model are presented in references [11] and [12] 
and will not be repeated here. The formulation presented in previous refer- 
ences is adapted to use the decohesion elements together with 2D plane strain 
elements. 

The constitutive model prevents interpenetration of the faces of the crack 
during closing, and a Fracture Mechanics-based criterion is used to predict 
crack propagation. The formulation of the damage model is summarized in 
Table 3, where and are the free energy per unit surface of the damaged 
and undamaged interface, respectively. % is the Kronecker delta, and d is 
a scalar damage variable. The parameter A is the norm of the displacement 
jump tensor (also called equivalent displacement jump norm), and it is used to 
compare different stages of the displacement jump state so that it is possible to 
define such concepts as ‘loading’, ‘unloading’ and ‘reloading’. The equivalent 
displacement jump is a non-negative and continuous function, defined as: 


A = \J (A,) 2 + (Aj) 2 (16) 

where (•) is the MacAuley bracket defined as ( x ) = \ (x + |x|). Ai is the 
displacement jump in mode I, i.e., normal to midplane, and A 2 is the dis- 
placement jump in mode II. 

[Table 3 about here] 

The evolution of damage is defined by G(-) , a suitable monotonic scalar func- 
tion, ranging from 0 to 1. /i is a damage consistency parameter used to define 


15 



loading-unloading conditions according to the Kuhn- Tucker relations. r x is the 
damage threshold for the current time, t, and r° denotes the initial damage 
threshold. The initial damage threshold is obtained from the formulation of 
the initial damage surface or initial damage criterion. A 0 is the onset displace- 
ment jump, and it is equal to the initial damage threshold r°. is the final 
displacement jump, and it is obtained from the formulation of the propagation 
surface or propagation criterion. Further details regarding the damage model 
used here can be found in references [11] -[12], 


4 Simulation of the double cantilever beam specimen 

The influence of mesh size, interface stiffness, interface strength, and the selec- 
tion of the parameters of the consitutive equation according to the proposed 
methodology were investigated by analyzing the Mode 1 delamination test 
(double cantilever beam, DCB). The DCB specimen was fabricated with a 
unidirectional T300/977-2 carbon- fiber-reinforced epoxy laminate. The spec- 
imen was 150-nnn-long, 20.0 nun-wide, with two 1.98-mm-thick arms, and it 
had an initial crack length of 55mm. The material properties are shown in 
Table 4 [32], 


[Table 4 about here] 


The FEM model is composed of two layers of four-node 2D plane strain ele- 
ments connected together with four-node decohesion elements. The decohesion 
elements were implemented using a user-written subroutine in the finite cle- 
ment code ABAQUS [33]. 
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Three sets of simulations were performed. First, a DCB test was simulated 
with different levels of mesh refinement using the material properties shown in 
Table 4 and interfacial stiffness of K=10 6 N/mm 3 . Then, equations (9) and (7) 
were used to calculate an adjusted interfacial strength and interface stiffness. 
Finally, a set of simulations with a constant mesh size using different interface 
stiffnesses were performed in order to investigate the influence of the stiffness 
on the calculated results. 

To study the effects of mesh refinement, several analyses were carried out for 
element sizes ranging between 0.125mm and 5mm. The corresponding load- 
displacement curves are shown in Figure 7. The results indicate that a mesh 
size of l e < 0.5mm is necessary to obtain converged solutions. The predictions 
made with coarser meshes overpredict significantly the experimental results. 
The length l cz of the cohesive zone for the material given in Table 4 is close 
to 1mm. Therefore, for a mesh size greater than 0.5mm, fewer than two el- 
ements would span the cohesive zone, which is not sufficient for an accurate 
representation of the fracture process [29]- [30]. 

[Figure 7 about here] 


4-1 Effect of interface strength 

To verify the effect of the interface strength on the computed results, simu- 
lations were performed by specifying the desired number of elements within 
the cohesive zone to be Nq = 5 and reducing the interface strength according 
to equation (11). The load-displacement curves obtained for several levels of 
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mesh refinement are shown in Figure 8. Accurate results are obtained for mesh 
sizes smaller than 2.5mm. 


[Figure 8 about here] 


[Figure 9 about here] 


A comparison of the load-displacement curves for the DCB specimen calcu- 
lated using the nominal interface strength and the strength obtained from 
equation (11) is shown in Figure 9. The specimen’s failure load obtained by 
keeping the maximum interfacial strength constant increases with the mesh 
size. Mesh dependence is especially strong for mesh sizes greater than 2mm. 
However, the failure loads predicted by modifying the interfacial strength ac- 
cording to equation (11) are nearly constant for element sizes smaller than 
3nnn. In addition, the global deformation and the crack tip position are also 
nearly independent from mesh refinement, as illustrated in Figure 10. 

[Figure 10 about here] 


J h 2 Effect of interface stiffness 

The DCB test was simulated with a mesh size of 2.5mm for various values 
of the interface stiffness in order to investigate the influence of the stiffness 
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on the predicted failure load. The results of the simulations are presented in 
Figure 11. 

The load-displacement response curves obtained from simulations using an 
interface stiffness greater than 10 4 N/mm 3 are virtually identical. However, 
smaller values of the interface stiffness have a strong influence on the load- 
displacement curves, since a stiff connection between the two neighboring lay- 
ers is not assured. Moreover, the number of iterations needed for the solu- 
tion when using an interface stiffness smaller than 10 4 N/nnn 3 is greater than 
the number of iterations needed for a range of the interface stiffness between 
10 6 * and 10 10 N/mm 3 . For values of the interface stiffness significantly greater 
than 10 10 N/mm 3 , the number of iterations needed for the solution increases 
(see Figure 12). The stiffness that results from equation (7) with a = 50 is 
5.55xl0 5 N/mm s , which is ideal for good convergence of the solution procedure. 


[Figure 11 about here] 


[Figure 12 about here] 


5 Simulation of free-edge delamination 

In order to demonstrate the mesh size effect and to further validate the accu- 

racy of the methodology proposed, simulations were performed of the fracture 
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of specimens without initial cracks. One of the causes of the initiation of de- 
lamination in specimens without initial cracks is the effect of free edge stress 
distributions. Delamination induced by the free-edge effect was investigated 
by many authors [34]- [36]. 


5.1 Problem statement 

Crossman et al. investigated the free edge effects in a [+25/-25/90] s graphite- 
epoxy laminate subjected to uniaxial strain [35]. The laminates were fabricated 
from Fiberite T300/934 prepreg tape using an autoclave. The cross-section 
of the specimen was 25.0-mm-wide and 0.792-mm-thick (see Figure 13). No 
initial cracks were induced in the laminate. The material properties are shown 
in Table 5 [34], [37]. 

Tensile tests were performed by applying a controlled displacement in the 
x-direction. Due to the stacking sequence, the through the thickness free- 
edge stress distribution produces high a zz in the 90° ply. Therefore, important 
tractions occur at the 90°/90° interface, leading to the onset of delamination 
at the free-edge. 


[Figure 13 about here] 


[Table 5 about here] 
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5.2 Numerical predictions 


A FEM model was developed using six layers of 4-node generalized plane strain 
elements [38]. The 90°/90° interface was modelled with 4-node decohesion 
elements presented in previous sections. The elastic properties of each layer 
were defined by means of the 3D stiffness matrix. The decohesion elements 
were implemented using a user-written subroutine in the finite element code 
ABAQUS [33]. Due to the symmetry of the laminate, only half of the cross- 
section was modelled. 

Two sets of simulations were performed. First, tensile tests were simulated 
with different levels of mesh refinement. Then, equations (9) and (7) were 
used to calculate the interface strength and the interface stiffness. 

Several analyses were carried out for mesh sizes ranging between 0.05nnn 
and 0.5mm. The load-displacement curves obtained are shown in Figure 14. 
An interfacial stiffness of K=10 6 N/mm 3 , an interface strength of 51.7MPa, 
a critical energy release rate of 0.175N/nnn [34], and the material properties 
shown in Table 5 were used. The length of the cohesive zone obtained with 
these properties, using equation (8), is 0.361mm. A mesh size smaller than 
0.1mm is needed in order to have three or more elements in the cohesive 
zone. The predictions obtained using different levels of mesh refinement are 
shown in Figure 14. Using meshes coarser than 0.1mm, the correct solution is 
significantly overpredicted. For a mesh size greater than 0.5mm the maximum 
applied load is greater than two times the maximum load obtained with meshes 
smaller than 0.1mm. 


[Figure 14 about here] 
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The ultimate tensile strength calculated using a mesh size of 0.1 mm is 404MPa. 
Delamination onset occurs when a point in the interface is not able to carry 
any traction, which changes the stress distribution. The stress maps in a plane 
perpendicular to the interface and normal to the load direction are represented 
in Figure 15 for three stages of deformation. At the first stage, before dclam- 
ination onset, the stresses near the free edge are tensile due to the mismatch 
effect of the Poisson ratio (black zones in Figure 15). At delamination on- 
set, a region of the free edge is unable to carry load and the stresses become 
compressive, as shown in Figure 15b). The predicted onset of delamination is 
e xx — 0.6%, and the corresponding stress is 392MPa. The experimental results 
obtained by Crossman and Wang [35] reported a delamination onset tensile 
stress of 409MPa and an ultimate tensile strength of 459MPa. The predicted 
and experimental onset and failure stresses and strains are summarized in 
Table 6. Unstable delamination propagation and the corresponding structural 
collapse occurs at a strain of 0.6183%. 

[Figure 15 about here] 

[Table 6 about here] 

To verify the effect of interface strength on the predicted results, simulations 
were performed by specifying the desired number of elements spanning the 
cohesive zone to be N 0 = 5 and reducing the interface strength according 
to equation (11). The load-displacement curves obtained for several levels of 
mesh refinement are shown in Figure 16. A comparison of the maximum load 
obtained using the nominal interface strength and the strength obtained from 
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equation (11) is shown in Figure 17. Accurate predictions are obtained using 
the modified interface strength instead of the nominal interface strength. 

[Figure 16 about here] 

[Figure 17 about here] 

A similar study of mesh size effect was repeated by the desired number of 
elements spanning the cohesive zone to be No = 10. Although the results 
presented in Figures 18 and 19 are more accurate than those for No = 5, the 
improvement is insignificant. 

[Figure 18 about here] 

[Figure 19 about here] 


6 Concluding remarks 

An engineering solution for the simulation of delamination using coarse meshes 
was developed. Two new guidelines for the selection of the parameters for the 
constitutive equation used for the simulation of delamination were presented. 
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First, a new equation for the selection of the interface stiffness parameter K 
was derived. The new equation is preferable to previous guidelines because 
it results from mechanical considerations rather than from experience. The 
approach provides an adequate stiffness to ensure a sufficiently stiff connec- 
tion between two neighboring layers, while avoiding the possibility of spurious 
oscillations in the solution caused by overly stiff connections. 

Finally, an expression to adjust the maximum interfacial strength used in the 
computations with coarse meshes was presented. It was shown that a minimum 
number of elements within the cohesive zone is necessary for accurate simula- 
tions. By reducing the maximum interfacial strength, the cohesive zone length 
is enlarged and more elements span the cohesive zone. The results obtained 
by reducing the maximum interfacial strength show that accurate results can 
be obtained with a mesh ten times coarser than by using the nominal inter- 
face strength. The drawback in using a reduced interfacial strength value is 
that the stress concentrations in the bulk material near the crack tip are less 
accurate. 
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Figures 



Fig. 1. Body 0 crossed by a material discontinuity T<i in the undeformed configu- 
ration. 
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Fig. 5. Length of the cohesive zone. 
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Fig. 6. Constitutive equations under Mode I and Mode II loading. 
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Fig. 7. Load- displacement curves using the nominal interace strength (t° =60MPa) 
for a DCB test with different mesh sizes. 
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Fig. 8. Load-displacement curves obtained for a DCB test with different mesh sizes 
with the interface strength modified to keep N e >5. 
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Fig. 9. Maximum load obtained in a DCB test for two cases: a) with constant inter- 
facial strength, b) with interfacial strength calculated according to Eq. (11). 
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Fig. 10. Comparison of crack tip position for two different levels of mesh refinement. 
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Fig. 11. Influence of the interface stiffness on the load- displacement curves. 
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Fig. 12. Influence 


of the interface stiffness on the number of iterations. 
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Fig. 14. Load- displacement curves obtained for a free-edge test with different mesh 
sizes using nominal interface strength properties. 
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Fig. 15. Evolution of stresses during delamination onset and propagation. 
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Fig. 16. Load- displacement curves obtained for a free-edge test with different mesh 
sizes and with the interface strength adjusted for N e =5. 


48 



Fig. 17. Maximum load obtained in a free- edge test for two cases: a) with constant 
interfacial strength, b) with interfacial strength calculated according to Eq.(ll) with 
N e =5. 
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Fig. 18. Load- displacement curves obtained for a free-edge test with different mesh 
sizes and with the interface strength adjusted for N e =10. 
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Fig. 19. Maximum load obtained in a free- edge test for two cases: a) with constant 
interfacial strength, b) with interfacial strength calculated according to Eq.(ll) with 
N e =1 0. 
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Tables 


Table 1 


Interface stiffness K proposed by different authors ( N/mm 3 ) and those calculated 
from equation (7). 


t(mm ) 0.125 

1 

2 

3 

5 

Equation (7) 4.43xl0 6 

5.5xl0 5 

2.75x10 s 

1.83xl0 5 

1.1x10 s 

Zou et al. [24] 

4.5xl0 5 

< K < 

4.5xl0 8 


Canranho and Davila [9] 10 6 

10 6 

10 6 

10 6 

10 6 
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Table 2 


Length of the cohesive zone and equivalent value of the parameter M. 



lez 

M 

Hui [26] 

2 p G c 

377 (r 0 ) 2 

0.21 

Irwin [25] 

1 E Gc 
* (r 0 ) 2 

0.31 

Dugdale [13], Barenblatt [14] 

TT 1 7* Gc 

8 ^( 7 ^ 

0.4 

Rice [27], Falk [21] 

97T if Gc 
32^(70)^ 

0.88 

Hillerborg [15] 

T? G c 

1 
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Table 3 


Definition of the constitutive model. 


Free Energy 


(A, d) = (1 - d) (Aj) - d i/j° (6u (-Ai)) 


Constitutive equation 


Ti = g = (i - d) SijKAj - d SijKSu (-Ai) 


Displacement jump norm A 


^/{A 1 ) 2 + (A 2 ) 2 


Damage criterion 


F (A 1 , ?’ 1 ) := G (A 1 ) - G (r*) < 0 Vt > 0 


G(A) 


a^(a~a°) 

A(A/-A») 


Evolution law 


_ :.dF{ A,r) _ - 9G(A) 


d = /i 


d\ 


= I*- 


dX 


r = jjL 


Load/unload conditions /i > 0 ; F (A 1 ,^) < 0 ; jiF (A^r 1 ) = 0 


r l = max { r° , max s A s } 0 < s < t 
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Table 4 


Mechanical and interface material properties of T300/977-2 [32]. 


Eu 

E22 — E33 

G 12 = G 13 

G23 

150.0GPa 

ll.OGPa 

G.OGPa 

3.7GPa 

3; 

to 

II 

3 

CO 

v 23 

Gic 

T° 

r 3 

0.25 

0.45 

0.352N/mm 

GOMPa 
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Table 5 


Mechanical properties of T300/934 Graphite Epoxy [34], [37]. 


E u 

CO 

CO 

II 

CN 

Gl2 = Gi3 

G 23 PL2 — ^13 

^23 

140.0GPa 

ll.OGPa 

5.5GPa 

3.61GPa 0.29 

0.52 
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Table 6 


Summary of predicted and experimental values of tension tests on T300/934 lami- 
nates, 



Predicted 

Experimental 

Difference 

Strain at delanrination onset (%) 

0.6 

0.59 

2% 

Stress at delanrination onset (MPa) 

392 

404 

-3% 

Ultimate strain (%) 

0.6183 

0.66 

-6% 

Ultimate strength (MPa) 

409 

459 

-10% 
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